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SUMMARY 

Calculations of steady and unsteady, transonic, turbulent channel flows 
with a time-accurate, lower-upper (L-U) factorization scheme are presented. 

The L-U factorization scheme is formally second-order accurate in time and 
space, and it is an extension of the steady state-flow solver (RPLUS) used 
extensively to solve compressible flows. A time-discretization method and the 
implementation of a consistent boundary condition specific to the L-U factori- 
zation scheme are also presented. The turbulence is described by the Baldwin- 
Lomax algebraic turbulence model. The present L-U scheme yields stable 
numerical results with the use of much smaller artificial dissipations than 
those used in the previous steady-flow solver for steady and unsteady channel 
flows. The capability to solve time-dependent flows is shown by solving very 
weakly excited and strongly excited, forced oscillatory, channel flows. 


INTRODUCTION 

In the lower-upper (L-U) factorization scheme (ref. 1), the governing 
flow equations are factored into forward (lower) and backward (upper) differ- 
ence operators, and the factored flow equations are solved by an implicit 
finite volume method. The L-U scheme has a few advantages over other numeri- 
cal methods for solving compressible flows. In the L-U scheme, only two 
factored equations need to be solved for any number of physical dimensions, 
whereas in the Alternating-Direction-Implicit (ADI) schemes, three factored 
equations need to be solved for three-dimensional flows. Thus, the number of 
numerical operations and the required memory size are smaller in the L-U 
schemes. Due to the use of an implicit time-stepping, which is usually free 
of the Courant-Friedrichs-Lewy (CFL) stability constraint condition, a 
relatively large time-step can be used in the L-U schemes. Even though the 
method is implicit, the costly matrix inversion of a large discrete system of 
equations is minimized, as shown in the following section. Due to the 
efficiency of the L-U scheme to solve large problems, the method has been 
rapidly extended to solve chemically reacting flows (ref. 2) and hypersonic 
flows (ref. 3). The method has also been extended to solve incompressible 
flows (ref. 4). 

Many fluid flows exhibit natural oscillations which may or may not be 
desirable from the design point of view. In certain cases, fluid flows may go 
through forced oscillations to achieve certain objectives. The start-up and 
shut-down of a fluid machinery also impart unsteadiness to the fluid flow, and 
critical fluid loading may occur during such processes. Therefore, unsteady- 
flow solvers are more desirable than steady-flow solvers. Furthermore, steady 
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flows can also be solved by using unsteady-flow solvers. In this report , the 
L-U scheme is extended to solve transient flows by incorporating the general- 
ized time-differencing scheme of Beam and Warming (refs. 5 and 6). To 
correctly resolve the transient flow phenomenon at and near the boundary , the 
boundary conditions are implemented in predictor and corrector steps in a 
mathematically consistent form with the time-accurate L-U factorization 
scheme. 


SYMBOLS 


A 

A 

B 

E,F 
E , F 

V V 

H 

I 

k 

(l,m) 

P 


Q 

R 


i 

X J 

As 


d 

K 

h 

A 


Jacobian matrix for £, A = {96/0Q} 

Jacobian matrix for F, B = {9 f/9q} 
convective transport vector 
diffusive transport vector 
height of channel throat , (H = 0.044 m) 
identity matrix 

effective thermal conductivity 

index for mesh 

pressure 

pressure at the exit boundary 
total pressure 

vector of flow variables, Q = {p , pu , pv , pe} 1 
pressure ratio, R p = P e /p t 
velocity component, u^ = {u,v> 

Cartesian coordinates, x = {x,y} 


intermediate solution 
phase angle 
molecular viscosity 
effective viscosity 
turbulent viscosity 
curvilinear coordinates, £ = 

order of magnitude 
Superscripts : 
n iteration count 

vector or matrix defined on curvilinear coordinates 
Subscripts: 

i,j index for spatial coordinates, i = {1,2,3} and j = {1,2,3} 
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L-U FACTORIZATION SCHEME 


The time-dependent, compressible, turbulent flow equation is given as 

9q _ _ 9e _ 9f 9e v 0f v (1 

(jt cJx Uy 9x 9y 

where 

Q = {p, pu, pv, pe} T 
E - {/>u, pu 2 + p, puv, u (pe + p)> T 
F = {pv, puv, pv 2 + p, v{pe + p)} T 

E v - < 0 ' T xx' T xy' UT xx + VT xy " ^ 

F v =* {0, T yx , T yy , u + vr yy - q y > T 

In the above equations, T = 2fi c + Ae 5 ^he stress tensor, 

e., - (Ou ,/0x. + OuVOx,) is the strain rate tensor, U = U + U is the effec- 

1 j X J • J 1 . 0 t 

tive viscosity, A = -2/3^ e is the second viscosity for compressibility, 
q i * -k0T/0x i is the heat flux due to temperature gradient, and the sub- 
scripts i and j denote each coordinate direction. For flows with arbi- 
trary geometries, it is more convenient to solve the flow equation defined on 
curvilinear coordinates with the use of body-fitted grids than to solve the 
flow equation defined on cartesian coordinates. The coordinate transformed 
flow equation can be written as 


9q ^ _ 9£ 9f 9fi v 9e v (2 ) 

9 r Hrj "5^ ~brj 

where Q = Q/J, fi - <f x E + £ y F + £ t Q)/J, F = (J? X E + 7/ y F + 7} t Q)/J, - (£ x E y + £ y F y )/J, 

F y = (ty x E v + T/yF v )/ J, J is the grid transformation Jacobian, and £ , f}^ r 

and 7jy are grid transformation metrics. A general form time-discretization 

for 0Q/0T 1 b given (ref. 5) as 


Aq" _ 9 i 9Aq" + l 9g n ^ e z Aq"’ 1 
"Sr l + W 2 9 r 1 + 0 2 9r 1 + W 2 br 


+ *\ 


- i - 4 


T + (AT)' 


(3) 


where the last term represents the truncation error. For the implicit, three- 
point, backward differencing case used in the present study, 8=1,9= 1/2, 

the temporal truncation error is d(&T) 3 , Aq r =* Q n+1 - Q n , and equation (3) is 
an alternative form of 
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3<AQ) n - (Aq)"- 1 
2At 


AA ->A n+1 „A n An-1 

OQ _ 3Q - 4Q + q 

2St 

Solving equation (2) for 9 qV9t and 5AQ n /9r and substituting the results 
Into equation (3) yields 


Aq" + 


e^r 


l + 9, 


0Ae n 9 Af" 
+ 


l 9 £ 9 7 


= R 


(4) 


where 


f 9g n 0F n ' 

(1 + 0 X )A T 

0g" 9 f" 

V V 

{ " 5T 

* l.». 

+ "Wj 


0 x Ar 


l + 6, 


gg"' 1 9F n_1 

r*r + - ^tj 


* 2 
i + 


Aq" 


In deriving the residual function (R) , Afi^ and Af" are approximated as 

^ n-1 [ A n-1 ] 

and AF v F v - F v J. In the ADI schemes (refs. 5 and 6) 

using the same time-discretization method (ref. 5), the time-lag approximation 
is used only for the cross-derivatives of the viscous transport terms. In the 
present case f all the viscous transport terms are lagged by one time step to 
facilitate the L-U factorization. The errors involved in these approximations 
are ft(Ar) , and hence the formal accuracy of the scheme remains intact 
(ref. 5). The convective transport terms in equation (4) are nonlinear 
functions of the flow variables. To solve equation (4) by a numerical method, 
the nonlinear terms are linearized as 


Afi 


n-lf 


s n - iH 

V V ] 


Ae" = aAq" + O(Ar) 2 . 
Af" = bAq” + ft(Ar) 2 


(5) 


where A = 0g n /9Q and B = 0F n /0Q are the Jacobian matrices. Substituting 
equation (5) into equation (4) yields 


Aq + 


0jA r 
1 + 9 ., 


9a 9b 

M + 


Aq = r 


( 6 ) 


where the superscript n has been deleted for notational convenience. For 
the L-U factorization of equation (6), the Jacobian matrices are decomposed as 
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(7) 


A 

A 


« A + A 

' + - 
* A + A 


where A + = 1/2 (A + rj), A' - 1/2 (A - r A I), B + = 1/2 (B + r B I), 

B“ = 1/2 (B - r fi I), r A > maxfji/J) and r B > max(|i/ B |), and and v are 

the eigen values of the Jacobian matrices A and B, respectively. The eigen 
values of^ A and B + matrices are positive values, and the eigen values of 
A and B matrices are negative values. Inserting equation (7) into 
equation (6) yields 


9 1 Ar 

( 

a~ 4+ 

b~, + 

9 + ,_ a",-) 

1 + 7T 

r A 

+ -5-® 

orj 

+ A +~ B 

9 £ ty J 


A Q = R 


( 8 ) 


where 9(A + + A )/9£ has been approximated as 0“A + /3£ + 9 + A’/9£, 0~A + /9£ 

- A (< ,m) - A + ( $ - l,m) is the backward difference operator and 0 + A~/3£ 

= A (< + l,m) - A (f,m) is the forward difference operator. The same approxi- 
mation is also applied to 9(B + + B )/9 7j . After some rearrangement, equa- 
tion (8) becomes 


0 X A T, 

1 + 7 r \^ + (^ m ) - A ( $ , m) -A + (f,m) + A (t +l,m) 

1 + & 2 


+ B ($ ,m) - B ($ ,m) - B + ($ ,m - 1) + B~(< ,m + 1) 


Aq = R 


or 


0 x At 

1 * TT 5 ; < ' : » * r,) 


9 Mr 


V x £17 r -I 

I - Jr A+ < { " 1 ' m > + - 1 ) 

1 + 1 


6 X Ar 

1 + ^ a- (r A + r B> 

1 + A B 


9 jAT r 1 

I + + 1 ' m ) + + 1)J 

1 + O'- 


Aq 


e^r 

“ni"-’ r «> 


R + 


e^r 


i + 9, 


V 


(9) 


<» - l,m) + B (f ,m - 1) 


x k~(( + l,m) + B (J ,m + 1 )]Aq 
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where A + - A~ * r^I and B + - B" = r B I have been made use of in deriving 
equation (9)* The last term in equation (9) represents the factorization 
error. The factorization error shares the same order of accuracy as those of 
various ADI factorization schemes (refs, 5 and 6), The L-U factored flow 
equation (eq. (9)) is solved in two steps: 


e^r 

1 + ttt 2 (r * + rB> 


r + 

1 ~ zr A 

nr 1 


l,m) + B (f ,m - 1) 


As 


:io) 


e x At 

1 + < r A + r B> 


i + 9, 


and 


[ M T 

9, At r _ i 

— T t lr * * V | 

I + [A ~ ( ? + 1 ,m) + B (*,m + l)j 

1 + 9 2 


Aq 


As 


(ID 


where As is the vector of intermediate solution. In the L-U factorization 
scheme, only two factored operators need to be solved for any number of 
physical dimensions; whereas in ADI schemes, three factored operators need to 
be solved for three-dimensional flows. Thus, less memory and computer time 
are required in the L-U schemes than in the ADI schemes. For flow equations 
without source terms (e.g., without chemical reactions), equations (10) 
and (11) are solved by forward and backward sweeps only; and the costly matrix 
inversion of a large discrete system of equations is unnecessary. 


CONSISTENT BOUNDARY CONDITION 

Usually, the boundary conditions in flow problems are given in terms of 
velocity, pressure and temperature; whereas the flow variables solved for in 
the present L-U scheme are the density, velocity, and internal energy. Thus, 
the boundary conditions given in terms of velocity, pressure, and temperature 
need to be recast in terms of density, velocity, and internal energy. Also 
note that the boundary conditions to be supplied in the lower sweep (see 
eq. (10)) need to be prescribed in terms of the intermediate solution on the 
boundary grid points. To correctly resolve time-dependent boundary condi- 
tions, the boundary conditions are implemented in predictor and corrector 
steps in the present L-U scheme. 


Predictor Step for Lower-Sweep 

The grid layout near the lower-sweep boundary is shown in figure 1. The 
intermediate solution on the west boundary or on the south boundary is 
obtained by using equation (11) and the three incremental solutions at the 
grid points AQ({,m), Aq($ + l,m) f and A<2(f,m + 1) shown in figure 1(c). The 
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incremental solutions at the grid points (f + l,m) and (J,m + 1) are predicted 
as 


Aq<* 
Aq (t 


+ l,m) = Q n ( f + l,m) - Q n ” 1 ( ( + l,m) 
+ l,m) =Q n (t + l,m) - Q n_1 ( t + l,m) 


(12) 


The incremental solution at the grid point ({ ,m) , or the grid point b in 
figure 1(d), is obtained from the boundary conditions and the flow variables 
of the most current time level. For example, by using the dependent 
variables, the temperature can be expanded as 


□ 5t . 5t 

T = T + _A/> + 

dp r 9(/?u) 


A(/>u) 


0T 

9(/>v) 


A (pv) 


9t 

9(/>e) 


A (pe) 


(13) 


For fixed temperature on the wall, the above equation can be rewritten as 


0T 0T 9t 
9/> ' 9(/>u) ' 9 (yOv) ' 


9t 

d(pe) 


[A p, A(/>u), A(/?v), A(/>e)] T 


T° - 


(14) 


where T° is the prescribed temperature which may or may not depend on time. 
For prescribed temperature gradient on the boundary, for the second-order one- 
sided difference approximation case, equation (13) can be rewritten as 


0T 

9n 



+ 


T a> 



or 


3 
2 

4 
2 

1 


9t 

V 

0T 

V 

[ 0T 


+ — 
2 


0 / 


9t 

0T 

0T 

9(/ju) ' 

9 </>v) ' 

9 (pe) 

9t 

0T 

0T 

9(/?u) ' 

9(/!v) ' 

9 (p&) 

0T 

0T 

0T 

9(/>u) ' 

9 (pv ) ' 

9(/?e) 


[A p, A(/>u), A (/>v), A(/>e)£ 
[A P, A(/>u) , A (/>v), A (/>e)£ 


[A/I, A(/)U), A(/Jv) , A(/>e)£ 


= T 


3 n 4 n In 

_T - _T - _T 
2 2 c 2 


(15) 


where the first, second, and third terms of equation (15) are evaluated at 
grid points b, c, and a, respectively (see fig. 1(d)). Derivation of the 
discrete boundary conditions for velocity and pressure follows the same 
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procedure as that for temperature. The discrete boundary condition can be 
written as 


§ a A 2 a + & b^b + ? (16) 

where £ = 6(T,u, v ,p) /d { p , pu f pv , pe) is the Jacobian matrix. The incremental 
solution at the boundary grid point, Aq , is obtained as 

b 

AQ b - e’V - § a A$ a + § c Aq c ) = Aq ( ? , m ) (17) 


Substituting equations (12) and 17) into equation (11) yields the boundary 
condition to be used in equation (10). 


Predictor Step for Upper-Sweep and Corrector Step 

For the upper sweep, the boundary condition on the east and north 
boundaries are obtained by applying equation (17) for each grid point on the 
boundary. For each time level, after lower and upper sweeps, the flow vari- 
ables on the boundary are corrected by using equation (17). 


NUMERICAL RESULTS AND DISCUSSION 

The geometry of the two-dimensional channel flows considered in the 
present study is shown in figure 2. Measured data for steady flow cases can 
be found in references 7 to 9 . For R (= p /p ) = 0.82, a mild normal shock 
appears in the downstream region of the throat. The flow field remains almost 
steady, and flow separation is not observed (ref. 7). For R = 0.72, a 
stronger normal shock than that for R =0.82 appears in the P farther down- 
stream region of the throat, and shock-induced separation occurs on the bottom 
and top walls (ref. 7). The flow for R = 0.72 exhibits a mild self-excited 
oscillatory motion. In the forced oscillation case, a triangular rotor, 
located at x=13H downstream of the throat, is turned at 100 revolutions per 
second. Thus, the flow is subjected to a variable exit pressure with the 
frequency of 300 Hz (refs. 7 and 8). The root-mean-square value of the first 
harmonic oscillatory pressure at x = 8.65H is approximately 0.57 percent of 
the local static pressure (refs. 7 and 8). However, the exact time-dependent 
exit pressure is not given in references 7 and 8. Hence, in the present 
study, the transient flow calculations are made by using harmonically varying 
exit pressure, as is done in references 9 and 10. 

In numerical calculations, for both steady and unsteady flows, the inlet 
boundary is located at 4H upstream of the nozzle throat, and the exit boundary 
is located at 12JT downstream of the throat. The entire flow domain is 
discretized by 113x81 grid points with the concentration of grid points in the 
near-wall and throat regions. The mesh for the contoured nozzle is shown in 
figure 3. The first grid point from the wall is located at y + » 1 (i.e., 
y/H = 0.651xl0~ 3 at the throat). The grid size in the transverse direction is 
increased by a factor of approximately 1.2. The mean velocities, static 
pressure, and temperature are prescribed at the inlet boundary. The static 
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pressure (time-dependent static pressure for unsteady flow) is prescribed at 
the exit boundary, and the velocities and the temperature are obtained by 
extrapolation. On the solid walls, no slip boundary condition is used for 
velocities, and vanishing gradient boundary condition is used for the pressure 
and temperature. The discrete system of equations is obtained by central 
differencing, and the numerical instability is suppressed by incorporating an 
artificial dissipation (ref. 11). The second- and fourth-order artificial 
dissipation coefficients used are 1/16 and 1/1280, respectively. These arti- 
ficial dissipation coefficients are much smaller than those used in the L-U 
factorization scheme for steady flows (refs. 2 and 3) and those used in the 
ADI schemes (refs. 6 and 11). For steady flows, the present method yields 
stable numerical results for the CFL number up to approximately 2000. The 
maximum CFL number that can be used to obtain stable numerical results is 
approximately one order of magnitude greater than those used in other numeri- 
cal methods (refs. 1 and 6). The unsteady flows are solved using the CFL 
number of 20. The capability to obtain stable numerical results by using a 
very small artificial dissipation and a large time-step is attributed to the 
time-accurate L-U factorization scheme. For steady flows, the converged 
solutions are obtained in approximately 2000 time steps. 

The turbulence is described by the Baldwin-Lomax turbulence model 
(ref. 12). It is known that the algebraic turbulence model may not yield 
accurate numerical results for largely separated flows (ref. 13). The present 
numerical calculations also reproduce the generally known shortcomings of the 
turbulence model. These shortcomings are discussed in more detail later in 
this section. 

The calculated pressure and Mach number contours for R =0.82 are shown 
in figure 4, where the incremental pressure and the incremental Mach number 
between the contour lines are constant. In figure 5, the pressure distribu- 
tions on the bottom and top walls are compared with the measured data of 
reference 7 and the numerical results of reference 9. The present numerical 
results are in good agreement with the measured data in most of the flow 
region except near x » 4.5H on the top wall. The lower pressure distribution 
in this region is caused by a small recirculation bubble located on the top 
wall at x * 4.5H. The recirculation bubble was not observed in the experi- 
ment (ref. 7). In the numerical calculation, the reversed flow region is 
caused by the algebraic turbulence model, which cannot accurately predict the 
rapidly growing turbulence intensity along the curved surface. In 
reference 9, the turbulence was described by the k-(t/ turbulence model, which 
did not show the reversed flow region. 

The calculated mean velocity profiles are shown in figure 6. The 
calculated mean velocity profile at x = 1.73H is in good agreement with the 
measured data and the other numerical result (ref. 9). Near and in the 
downstream region of the spurious separation bubble, the present numerical 
results compare less favorably with the measured data than do the other 
numerical results . 

The calculated pressure and Mach number contours for R =0.72 are shown 
in figure 7. Again, the incremental pressure and the incremental Mach number 
between the contour lines are constant. It can be seen in the figure that the 
shock is moved further in the downstream direction and that the Mach number 
decreases more rapidly across the shock than the previous case (R =0.82) due 


9 



to the stronger shock. The calculated pressure distributions on the bottom 
and top walls are shown in figure 8. The calculated flow field contains a 
large shock-separated recirculation zone on the top wall and a small shock- 
separated recirculation zone on the bottom wall. The relatively large 
discrepancy between the calculated and the measured wall pressure distribu- 
tions is caused by the large internal blockage formed by the recirculation 
bubbles. In reference 10 , the exit boundary was located at 8.64H downstream 
of the nozzle throat. In their calculation using the original k-W turbu- 
lence model, the reversed flow region extended beyond the exit boundary and 
converged solution was not obtained. The converged solution was obtained by 
using a different turbulence model constant in reference 10. The use of a 
different turbulence model constant effectively reduces the size of the 
recirculation bubble and yields numerical results which are in good agreement 
with the measured data (ref. 10). The same turbulence model that was used in 
reference 10 has also been used in reference 9. In the present calculation 
using the Baldwin-Lomax turbulence model, a converged solution cannot be 
obtained iT the exit boundary is located at 8.64H downstream of the nozzle, 
since the reversed flow region extends beyond the exit boundary. The con- 
verged solution is obtained by locating the exit boundary at 12H downstream of 
the nozzle throat instead of optimizing the Baldwin-Lomax turbulence model. 

The calculated mean velocity profiles for R =0.72 are shown in fig- 
ure 9. The calculated mean velocity profiles show a large reversed flow 
region near the top wall boundary. Near and downstream of the separation 
region, the present numerical results compare less favorably with the measured 
data than do the numerical results of reference 9. The present numerical 
results, as well as those of reference 10, strongly demonstrate that an 
accurate numerical result can hardly be obtained without the use of a good 
turbulence model. However, the present study is limited to the development 
and verification of the time-accurate L-U factorization scheme, and optimiza- 
tion of the algebraic turbulence model is not attempted. 

The evolution of flow variables in time for mildly and strongly excited 
unsteady channel flows are shown in figures 10 and 11, respectively. The exit 
pressure is prescribed as 


P e = P 0 + P c sin(27Tft) 


where t represents time. For the mildly excited case, P =0.82, 

P c = 0.0082, and f = 300; and for the strongly excited case, p =0.77, 

P c = 0.05, and f = 300. In each case, the numerical results show that the 
fluid motion consists of the low-frequency natural oscillation (which is 
inherent to each flow system) and the forced oscillation. The oscillatory 
motion of the separation bubble for the strongly excited case is shown in 
figure 12. As the exit pressure is increased, the separation bubble becomes 
bigger and moves toward the throat; and as the exit pressure is decreased, the 
bubble becomes smaller and moves toward the exit boundary. Thus, the numeri- 
cal method reproduces the experimentally observed trend of the fluid motion. 
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SUMMARY OF RESULTS 


Numerical calculations of steady and forced oscillatory transonic 
turbulent flows in a two-dimensional channel are made by using a time-accurate 
L-U factorization scheme. The calculated velocity profiles and the pressure 
distributions on the walls for the steady flows are in good agreement with the 
measured data, considering the limited predictive capability of the algebraic 
turbulence model for complex turbulent flows. The present numerical results 
show a large shock-separated reversed flow region for the R =0.72 case. 

The large separation bubble is caused by the use of the alge&raic turbulence 
model which cannot adequately describe the turbulence field subjected to the 
streamline curvature and the normal shock. It is interesting to note that the 
k-W turbulence model (ref. 10) also yields a largely separated flow region for 
the steady flow for R =0.72. The objective of the present study is limited 
to the development and the verification of the time-accurate lower-upper 
factorization scheme, and optimization of the turbulence model is not 
attempted. 

The time-accurate lower-upper factorization scheme shows a few important 
numerical aspects. A large amount of artificial dissipation can impair the 
numerical results and hence alternative artificial dissipation models have 
also been proposed (refs. 14 and 15). However, the present numerical method 
yields stable numerical results by using a much smaller artificial dissipation 
than other numerical methods (refs. 2, 3, 6, and 11). The method also yields 
a stable numerical result with the use of a large CFL number. The maximum CFL 
number that can be used to obtain stable numerical results for the present 
method is at least one order of magnitude greater than those for the other 
numerical methods (refs. 1, 6, and 11). Thus, the method is competitive with 
steady flow solvers for steady flows. For time-dependent flows, the time-step 
size is limited only by the physical phenomenon of the fluid flow to be 
considered. 
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Figure 3. — Mesh in the vicinity of channel throat. 
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Figure 4. — Contour plot of pressure and Mach number for p e /p t = 0.82. 
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Figure 5. — Pressure distributions for p e /p t = 0.82. 
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Figure 6. — Steady-state velocity profiles for 
p e /p t = 0.82. Profile (a), x/H = 1.73; (b), x/H = 
4.04; (c), x/H = 6.34; (d) x/H = 8.36. 
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Figure 7. — Contour plot of pressure and Mach number for p 6 /p t - 0.72. 


15 




o Expt. (ref. 7) 

Hsieh (ref. 9) 

Present 



(a) Bottom wall. 



cl .70 

Q. 

¥ 65 

<0 

.70 
3 .65 


P (exit) 

p (x/H = 4.74, y/H = 0) 
p (x/H = 4.74, y/H = 1.43) 
p (x/H = 4.74, y/H = 0.72) 

u (x/H = 4.74, y/H =0.72) 


Figure 10. — Evolution of velocity and pressure 
in time for very weakly excited flow (p e = p 0 + 
p c sfn(2TrfT); where p 0 = 0.82, p c = 0.0082, 
and f = 300). 
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Figure 9. — Steady-state velocity profiles for 
Pe/Pt = 0.72. Profile (a), x/H = 2.89; (b). x/H : 
4.61 ; (c), x/H = 5.76; (d) x/H = 8.36. 
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Figure 11. — Evolution of velocity and pressure 
in time for strongly excited flow (p e = p 0 + 
p c sin(27rfT); where p 0 = 0.77, p c = 0.05, and 
f = 300). 
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Figure 12. — Oscillatory motion of separated bubble for strongly excited flow. 
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